CRITICAL PARAMETER VALUES AND RECONSTRUCTION PROPERTIES OF DISCRETE 
TOMOGRAPHY: APPLICATION TO EXPERIMENTAL FLUID DYNAMICS 

STEFANIA PETRA, CHRISTOPH SCHNORR, ANDREAS SCHRODER 

Abstract. We analyze representative ill-posed scenarios of tomographic PIV with a focus on conditions for 
unique volume reconstruction. Based on sparse random seedings of a region of interest with small particles, 
the corresponding systems of linear projection equations are probabilistically analyzed in order to determine 
(i) the ability of unique reconstruction in terms of the imaging geometry and the critical sparsity parameter, 
and (ii) sharpness of the transition to non-unique reconstruction with ghost particles when choosing the sparsity 
'qj ' parameter improperly. The sparsity parameter directly relates to the seeding density used for PIV in experimental 

fluids dynamics that is chosen empirically to date. Our results provide a basic mathematical characterization 
of the PIV volume reconstruction problem that is an essential prerequisite for any algorithm used to actually 
compute the reconstruction. Moreover, we connect the sparse volume function reconstruction problem from 
few tomographic projections to major developments in compressed sensing. 
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1. Introduction 



< 

Motivated by an application from fluid dynamics [9], we investigate conditions for an highly underdeter- 
S mined nonnegative system of linear equations to have a unique nonnegative solution, provided it is sparse. 

The sought solution is a sparse 3D image of particles immersed in a fluid known only from its projec- 
i-H tion. This projection represents the simultaneous 2D images captured by few camera sensors from different 

^ viewing directions, see Fig. 1. The reconstruction of the 3D image from the 2D images employs a standard 

algebraic image reconstruction model, which assumes that the image consists of an array of unknowns (cells, 
^ voxels), and sets up algebraic equations for the unknowns in terms of measured projection data. The latter 

are the pixel entries in the recorded 2D images that represent the integration of the original 3D light intensity 
distribution along the pixels line-of-sight. The number of cameras is limited to 3 to 6 cameras, typically 4. 
O As a consequence, the reconstruction problem becomes severely ill-posed. 

Thus, we consider a huge and severely underdetermined linear system 



Ax^h, AeM^^^, m<n, (1.1) 

with the following properties: a very sparse nonnegative measurement matrix A with constant small support 
of length i of all column vectors, 

A > 0, X > 0, supp(A.j) = ^ < m, Vj = 1, . . . , n. (1.2) 

and a nonnegative fc-sparse solution vector x. While i equals the number of cameras, k is equal related to 
the particle density (equal, in the present work, proportional, in practice). We also consider the discretization 
(or resolution) parameter d, and relate it to the number of discretization cells and number of measurements: 

m = £-0{d), n^O{d^), in 2D, (1.3) 
m = £-0{d^), n^O{d^), in 3D . (1.4) 



Key words and phrases, compressed sensing, underdetermined systems of linear equations, sparsity, large deviation, tail bound, 
algebraic reconstruction, TomoPIV. 
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Figure 1 . Typical camera arrangements: in circular configuration (right) or all in line (left). 



We will answer the following question: what is the maximal number of particles, depending on the image 
resolution parameter d, that can be localized uniquely? Formally, we want to relate the exact recovery of x 
from it's noiseless measurements b to the sparsity k and to the dimensions of m, n of the projection matrix 
A. Moreover, we investigate critical values of the sparsity parameter k such that most fc-sparse nonnegative 
solutions are the unique nonnegative solutions of (1.1) with high probability. 

1.1. Related Work. Research on compressed sensing [5, 3] focuses on properties of underdetermined linear 
systems that guarantee exact recovery of sparse or compressible signals x from measurements b. Donoho and 
Tanner [7, 8] have computed sharp reconstruction thresholds for random measurement matrices, such that 
for given a signal length n and numbers of measurements m, the maximal sparsity value k which guarantees 
perfect reconstruction can be determined explicitly. The authors derived their results by connecting it to 
problem from geometric probability that n points in general position in can be linearly separated [15]. 
This holds with probability Pr(n,m) = 1 for n/m < 1, and with Pr(n,m) ^ 1 if m ^ oo and 1 < 
n/m < 2, where 

m— 1 



(1.5) 



The authors show in [7, Thm. 1.10] that the probability of uniqueness of a fc-sparse nonnegative vector 
equals Pr(n — m, n — fc), provided A satisfies certain conditions which do not hold in our considered appli- 
cation. Likewise, by exploiting again Wendel's theorem, Mangasarian and Recht showed [11] that a binary 
solution is most likely unique if m/n > 1/2, provided that A comes from a centrosymmetric distribu- 
tion. Unfortunately, the underlying distribution A lacks symmetry with respect to the origin. However, we 
recently showed [13] for a three camera scenario that there are thresholds on sparsity (i.e. density of the par- 
ticles), below which exact recovery will succeed and above which it fails with high probability. This explicit 
thresholds depend on the number of measurements (recording pixel in the camera arrays). The current work 
investigates further geometries and focus on an average case analysis of conditions under which uniqueness 
of X can be expected with high probability. A corresponding tail bound implies a weak threshold effect and 
criterion for adequately choosing the value of the sparsity parameter k. 

1.2. Notation. |X| denotes the cardinality of a finite set X and [n] = {1, 2, . . . , n} for n G N. We will 
denote by ||x||o = xi ^ 0}| and = {x G M^: ||x||o < k} the set of fc-sparse vectors. The 
corresponding sets of non-negative vectors are denoted by W]^ and M^^, respectively. The support of a 
vector X E M^, supp(x) C [n], is the set of indices of non- vanishing components of x. 

For a finite set S, the set A/'(S') denotes the union of all neighbors of elements of S where the correspond- 
ing relation (graph) will be clear from the context. 
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A.^i denotes the z-th column vector of a matrix A. For given index sets /, J, matrix Ajj denotes the 
submatrix of A with rows and columns indexed by / and J, respectively. I^, denote the respective 
complement sets. Similarly, 6/ denotes a subvector of b. 

E[-] denotes the expectation operation applied to a random variable and Ft{A) the probability to observe 
an event A. 



2. Graph Related Properties of Tomographic Projection Matrices 

Recent trends in compressed sensing [2, 16] tend to replace random dense matrices by adjacency matrices 
of "high quality" expander graphs. Explicit constructions of such expanders exist, but are quite involved. 
However, random m x n binary matrices with nonreplicative columns that have [in\ entries equal to 1, 
perform numerically extremely well, even if £ is small, as shown in [2]. In [10, 12] it is shown that perturbing 
the elements of adjacency matrices of expander graphs with low expansion, can also improve performance. 

2.1. Preliminaries. For simplicity, we will restrict on situations were the intersection lengths of projection 
rays corresponding to each camera with each discretization cell are all equal. Thus, we can make the as- 
sumption that the entries of A are binary. It will be useful to denote the set of cells by C = [n] and the set of 
rays by i? = [m] . The incidence relation between cells and rays is then given by 

, , , I 1, if 7-th ray intersects i-th cell, 

{A)ij = { ' ^ ' (2.1) 

10, otherwise, 

for all z G [m], j G [n]. Thus, cells and rays correspond to columns and rows of A. 

This gives the equivalent representation in terms of a bipartite graph G = (C, R] E) with left and right 
vertices C and i?, and edges (c, r) ^ E iff {A)rc = 1. G has constant left-degree i equal to the number of 
projecting directions. 

For any non-negative measurement matrix A and the corresponding graph, the set 

X{S) = [m] : Aj > 0, j G S} 

contains all neighbors of S. The same notation applies to neighbors of subsets S C [m] of right nodes. 
Further, we will call any non-negative matrix adjacency matrix, based on the incidence relation of its non- 
zero entries. 

If A is the non-negative adjacency matrix of a bipartite graph with constant left degree i, the perturbed 
matrix A is computed by uniformly perturbing the non-zero entries Aij > to obtain Aij G [Aij —5, Aij +5] , 
and by normalizing subsequently all column vectors of A. In practice, such perturbation can be implemented 
by discretizing the image by radial basis functions of unequal size or and choose their locations on an 
irregular grid. 

The following class of graphs plays a key role in the present context and in the field of compressed sensing 
in general. 

Definition 2.1. A (z/, 5)-unbalanced expander is a bipartite simple graph G = (L, R; E) with constant left- 
degree i such that for any X C L with |X| < z/, the set of neighbors J\f{X) C R of X has at least size 

\M{X)\ >6£\X\. 

Recovery of a fc-sparse nonnegative solution via an (z/, 5) -unbalanced expander was derived in [14]. It 
employs the smallest expansion constant S with respect to other similar results in the literature. 

Theorem 2.1. Let A be the adjacency matrix of a {y^ S)-unbalanced expander and 1 > 6 > - Then for 
any k-sparse vector x* with k < (^i^^y the solution set {x : Ax — Ax*, x > 0} is a singleton. 
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Now let A denote the tomographic projection matrix, and consider a subset X C C of \X\ = k columns 
and a corresponding fc-sparse vector x. Then b = Ax has support J\f{x), and we may remove the subset of 
Af{Xy = {J\f{X)y rows from the linear system Ax = b corresponding to 6^ = 0, Vr G i?. Moreover, 
based on the observation J\f{X), we know that 

X C AA(AA(X)) and AA(AA(X)") n X = 0. (2.2) 

We continue by formalizing the system reduction just described. 

Definition 2.2. The reduced system corresponding to a given non-negative vector 6, 

^red^ — i^redi ^red ^ M^^^"^ ^ (2.3) 

results from A, 6 by choosing the subsets of rows and columns 

i?5 := supp(6), C5 := AA(i?5) \M{Rl) (2.4) 

with 

^red •= nred •= iC'fel- (2.5) 

Note that for a vector x and the bipartite graph induced by the measurement matrix A, we have the 
correspondence (cf. (2.2)) 

X = supp(x), i?5 = AA(X), Cfe = M{N{X)) \N{N{Xf). 

We further define 

5+ := {x: Ax = > 0} (2.6) 

and 

^red ^i^bC'b^ = > 0} . (2.7) 

The following proposition asserts that solving the reduced system (2.3) will always recover the support of 
the solution to the original system Ax = b 

Proposition 2.2. [13, Prop. 5.1] Let A G W^^^ and b G have nonnegative entries only, and let and 
^'red defined by (2.6) and (2.7), respectively. Then 

5+ = {x G M^: X(c,)c = and xc, G 5+ J. (2.8) 

Consequently, we can restrict the linear system Ax = 6 to the subset of columns J\f{J\f{X)) \ 
J\f(J\f(Xy) C C, and only consider properties of this reduced systems. 

2.2. Guaranteed Uniqueness. Uniqueness of x G is guaranteed if all k or less-sparse supported on 

supp(x) induce overdetermined reduced systems with mred/^red > Si > ^^^ i. 

Proposition 2.3. [13, Th. 3.4] Let A he the adjacency matrix of a bipartite graph such that for all random 
subsets X C C of\X\ < k left nodes, the set of neighbors N {X) ofX satisfies 

\M{X)\ > 5i\M{M{X)) \M{M{xy)\ with 5 > ^^^^^ (2.9) 

Then, for any Sk-sparse nonnegative vector x*, the solution set {x : Ax — Ax*, x > 0} is a singleton. 

For perturbed matrices uniqueness is guaranteed for square reduced systems, and thus less high sparsity 
values. 



WEAK RECOVERY THRESHOLDS FOR TOMOPIV 



5 




Figure 2. Sketch of a 3-cameras setup in 2D. Left: The hexagonal area discretized in 
3^^^ equally sized cells projected on 3 ID cameras. The resulting projection matrix 

A G {0, is underdetermined, with m = 3d and n = 3^^, where ^ + 1 is 

the number of cells on each hexagon edge. Middle: This geometry can be easily extended to 
3D by enhancing both cameras and volume by one dimension, thus representing scenarios of 
practical relevance when cameras are aligned on a line. Right: When considering a square 
area along with three projection directions (two orthogonal, one diagonal) one obtains a pro- 
jection matrix with analogous reconstruction properties. The projection matrix equals up to 
scaling the previous projection matrix corresponding to the hexagonal area if we remove the 
2 • ^^g^ cells in the marked corners along with incident rays. 



Proposition 2.4. [13, Th. 3.4] Let A be the adjacency matrix of a bipartite graph such that for all subsets 
X C C of\X\ < k left nodes, the set of neighbors N {X) of X satisfies 

\M{X)\ > 5l\M{M{X)) \M{M{Xf)\ with 5>j. (2.10) 

Then, for any k-sparse vector x*, there exists a perturbation A of A such that the solution set {x: Ax — 
Ax*, X > 0} is a singleton. 

Recovery via perturbed underdetermined reduced systems is possible and our numerical results from 
Section 5 suggest the following. 

Conjecture 2.5. Let A be the adjacency matrix of a bipartite graph such that for all subsets X C C of 
\X\ < k left nodes, the set of neighbors M {X) ofX satisfies 

W{X)\ > l±^\M{M{X))\M{M{xr)\ with 5 > (2.11) 

Then, for any j -sparse vector x*, there exists a perturbation A of A such that the solution set {x: Ax = 
Jlx*, X > 0} is a singleton. 

The consequences of Propositions 2.3, 2.4 and Conjecture 2.5 are investigated in the following sections 
3.2 and 4.2 by working out critical values of the sparsity parameter k for which the respective conditions are 
satisfied with high probability. 

3. 3 Cameras - Left Degree equals 3 

In this section, we analyze the imaging set-up depicted in Figure 2, left panel, which also represents typical 
3D scenarios encountered in practice with a coarse resolution only along the third coordinate, as shown by 
Figure 2, center panel. 
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3.1. Imaging Geometry. Cell centers Xc of hexagonal cells c G C that partition a region of interest, are 
given by lattice points corresponding to integer linear combinations of two vectors rf^ z = 1, 2, 

x,^hd^ + i2d\ = (n,i2)ex, (3.1) 

for the index set 

1 = {(z, j) : - (d - l)/2 < i, j < (d - l)/2, \i + j\ <{d- l)/2}, (3.2) 

with problem size d G N that we assume (in this section) to be an odd number for simplicity. The number 
\R\ of projections r G i? = i?i U i?2 U Rs, where i?^, z = 1, 2, 3, corresponds to the rays of direction i, is 

\R\ = 3\R.\ = 3d. (3.3) 

The number of cells incident with projection rays ranges over the interval 

{(d+l)/2,(d+l)/2 + l,...,d} (3.4) 

from the periphery towards the center. Thus, indexing with r each projection ray along any particular direc- 
tion i = 1, 2, 3, from one side of the hexagon across the center towards the opposite side, the numbers 
of cells incident with ray r is 

\r\ E {{d + l)/2, . . . , d, . . . , (d + l)/2}, r G i?^, i = 1, 2, 3. (3.5) 

The total number of cells is 

d-l 

|C|=5^|r| = 2 Yl j + d=-{3d^ + l), zG {1,2,3}, (3.6) 

reRi j={d^l)/2 



and 



^|r| = 3|C|. (3.7) 



reR 

Accordingly, the system of equations representing the imaging geometry depicted by Figure 2, left panel, 
has dimensions 

Ax^b, AG {0,1}I^I^I^I, 6gMI^I. (3.8) 

Note that \R\ <C \C\. For further reference, we define the quantities 

\r\ 

q = mm qr, 

and list some further relations and approximations for large d, 

_ 2(d+ 1) _ 2 
^~ 3^2 + 1 ~ 3d' 
\R\q^^2, I^l9d~4. (3.10b) 



Pr = 


1 - ^r, 




(3.9a) 


Qd = 


maxg^, 

reR 




(3.9b) 


Pd = 






(3.9c) 


Qd 


1 

~3p' 


4 

' 3d' 


(3.10a) 
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3.2. Dimensions of Reduced Systems. We estimate the expected dimensions (2.5) of the reduced system 
(2.3) based on uniformly selecting k cells at random locations. 

To each projection ray r ^ R,wq associate a binary random variable Xr taking the value = 1 if not 
any of the k cells is incident with ray r, and = otherwise. We call the event = 1 zero-measurement. 

We are interested in the random variable 



r^R 

that determines the number of projection rays not incident with any of the k cells, that is the number of zero 
measurements. We set 



Hence, Nr is the expected size of the support rrired = I supp(6) | of the measurement vector h. 

Remark 3.1. Note that random variables X^ are not independent because different projection rays may 
intersect. This dependency does not affect the expected value of X, but it does affect the deviation of 
observed values of X from its expected value - cf. Section 3.3. 

Remark 3.2. We do not assume in the derivation below that k different cells are selected. In fact, a single cell 
may be occupied by more than a single particle in practice, because real particles are very small relative to 
the discretization cells c. The imaging optics enlarges the appearance of particles, and the action of physical 
projection rays is adequately represented by linear superposition. 

Definition 3.1 (Sparsity Parameter). We refer to the number k introduced above as sparsity parameter. 
Thus, highly sparse scenarios correspond to low values k. 

Lemma 3.1. The expected number of zero measurements is 




(3.11) 



iV° :=E[X], 



Nr := \R\ - N\ 



(3.12) 




(3.13) 



r&R 



Proof. For k — 1, Xr has a Bernoulli distribution with 



E[Xr] = Vx[Xr = 1] = 1 




(3.14) 



For k independent trials, we have (cf. Remark 3.2) 

E.[Xr] = Pr[X^ = 1] 
By the linearity of expectations and (3.15), we obtain (3.13), 




(3.15) 




(3.16) 



□ 



We discuss few specific scenarios depending on the sparsity parameter k. 
No particles: For A; = 0, we obviously have 



< = J^l = |i?|. 



reR 

High sparsity: By (3.10), we have < Qr < Qd^ hence qr = 0{d~^). Thus, for large problem sizes 
d and small values of k. 
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By (3.7), J2reR = 3, hence 

N^^ \R\ -3k. (3.17) 

This approximation says that for sufficiently small values of k each randomly selected cell can be 
expected to create 3 independent measurements, which just reflects the fact that each cell is met by 
three projection rays. 

Less high sparsity: For increasing values of k higher-order terms cannot longer be ignored, due to the 
increasing number of projection rays meeting several particles. Taking the second-order term into 
account, we obtain in an analogous way 

k(k-l) 3 ^^-'^^ 

< 5^(1 - kqr + A_dq^q^) ^ \R\ _ 3fc + -k(k - l)q,, 

which is a fairly tight upper bound for values of k and N that are relevant to applications. 

We consider next the expected number of cells rired = IC'fel of cells supporting the set i?^ according to 
(2.4). We denote this expected number by 

7Vc:=E[|C5|], N^:=\C\-Nc^ (3.19) 

and by the expected size of the complement. 

Let i? = i?i U i?2 U i?3 denote the partition of all projection rays by the three directions. For each cell 
c, there are three unique rays ri{c) G i = 1, 2, 3, incident with c. Furthermore, for i ^ j and some ray 
Ti ^ Ri, let Rj{ri) denote the set of rays that intersect with r^. As before, \r\ denotes the number of cells 
covered by projection ray r ^ R. 

Proposition 3.2. For a given sparsity parameter k, the expected number of cells that can be recognized as 
empty based on the observations of random variables {Xr}r^Ry 

= N^(k) = 3iV^ - 37V^ + N^, (3.20a) 

T^Ri 



^c=Y. E (i-^^^^^p— ^) ,>'-««y^,ie{i,2,3}, i^j, (3.20c) 

g|-._ ELMc)i-2 y ^^^^^ 



Proof Each cell intersects with three projection rays ri(c), i = 1, 2, 3. Hence, given the rays corresponding 
to zero measurements, each cell that can be recognized as empty if either one, two or three rays from the set 
{n(c)}i=i,2,3 belong to this set. 

We therefore determine separately the expected number of removable cells (i) due to individual rays cor- 
responding to zero measurements, (ii) due to all pairs of rays that intersect and correspond to zero measure- 
ments, and (iii) due to all triples of rays that intersect and correspond to zero measurements. The estimate 
(3.20a) then results from the inclusion-exclusion principle that combines these numbers so as to avoid over- 
counting, to obtain the desired estimate corresponding to the union of these events. 
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Consider each projection ray r ^ Ri for any fixed direction i = 1,2,3. Because these rays do not intersect, 
the expected number of cells that can be removed based on the observation {Xr}reR^ is 

Nh^E[Y,Xr\r\\ = J]p^|r|, (3.21) 

reRi reRi 

by the linearity of expectations and (3.15). Due to the symmetry of the setup, this number is the same for 
each direction z = 1, 2, 3. Hence we multiply by 3 in (3.20a). 

Consider next pairs of directions z,jG{l,2,3}, i^^j. For i fixed, the expected number of empty cells 
based on a zero measurement corresponding to some ray G Ri and all rays Vj G Rj{ri) intersecting with 

is 

Nc-^[Y. E ^n^r\ (3.22) 

Ti^Ri rjeRjiri) 

Thelinearity of expectations and E[X^.X^^ ] = Pr [(X^. = 1)A(X^^ = 1)] gives (3.20c). Due to symmetry, 
we have to multiply by 3 in (3.20a). 

Finally, the expected number of empty cells that correspond to observed zero measurements along all 
three projection directions, is 

3 

(3.23) 

'c^Ci=l ^ 

which equals (3.20d). □ 



An inmiediate consequence of Lemma 3.1 and Prop. 3.2 is 

Corollary 3.3. For a given value of the sparsity parameter k, the expected dimensions of the reduced system 
(2.3) are 

rrired = Nr- N^, Tired = Nc - X^, (3.24) 
with 7V^, given by (3.13) and (3.20). 

3.3. A Tail Bound. We are interested in how sharply the random number X of zero measurements peaks 
around its expected value = E[X] given by (3.13). 

Because the random variables Xr^ r e R, are not independent due to the intersection of projection rays, 
we apply the following classical inequality for bounding the deviation of a random variable from its expected 
value based on martingales, that is on sequences of random variables (X^) defined on a finite probability 
space (fi, J^, /x) satisfying 

E[Xi+i| j;] = Xi, for all i > 1, (3.25) 
where denotes an increasing sequence of cr-fields in T with Xi being -measurable. 

Theorem 3.4 (Azuma's Inequality [1, 4]). Let (X^)^=o,i,2,... be a sequence of random variables such that 
for each i, 

\X,-X,_i\<c,. (3.26) 

Then, for all j > and any 5 > 0, 

Pr {\Xj - Xol > (5) < 2 exp ( - ) . (3.27) 



Let C 2^, i = 0, 1, 2, ... , denote the cr-field generated by the collection of subsets of R that corre- 
spond to all possible events after having observed i randomly selected cells. We set = {0, R}- Because 
observing cell z + 1 just further partitions the current state based on the previously observed i cells by pos- 
sibly removing some ray (or rays) from the set of zero measurements, we have a nested sequence (filtration) 
^0 ^ ^1 ^ • • • ^ ^/c of the set 2^ of all subsets of R. 
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Based on this, for a fixed value of the sparsity parameter k, we define the sequence of random variables 

= E[X|j;], z = 0,l,...,A:, (3.28) 

where Y^, z = 0, 1, . . . , A: — 1, are the random variable specifying the expected number of zero measurements 
after having observed k randomly selected cells, conditioned on the subset of events determined by the 
observation of i randomly selected cells. Consequently, Yq = E[X] = due to the absence of any 
information, and = X is just the observed number of zero measurements. The sequence (li)i=o,...,fc is a 
martingale by construction satisfying E[y^+i| J'i] = Yi, that is condition (3.25). 

Proposition 3.5. Let — E[X] be the expected number of zero measurements for a given sparsity param- 
eter k, given by (3.13). Then, for any 5 > 0, 

^2 



Pr [|X - Nl\ >5] < 2exp f ] 

u R\ - \ - 18(1 -pf) J 



(3.29) 



Proof Let Ri_i C R denote the subset of rays with zero measurements after the random selection of 
i — 1 < k cells. For the remaining k — {i — 1) trials, the probability that not any cell incident with some ray 
r G Ri_i will be selected, is 

p^-(^-i) ^EfX^lJ-.-i], (3.30) 
with pr given by (3.14). Consequently, by the linearity of expectations, the expectation Yi-i of zero mea- 
surements, given the number of zero measurements after the selection of z — 1 cells, is 

Fi-i = = Yl Pr~^'~^^- (3-31) 

Now suppose we observe the random selection of the i-th cell. We distinguish two possible cases. 

(1) Cell i is not incident with any ray r e R^-i- Then the number of zero measurements remains the 
same, and 

Yi- P^r~'- (3-32) 

Furthermore, 

<Pd-'T.1r = ^Pd-'- 

(2) Cell i is incident with one, two or three rays contained in R^_i. Let R^ denote the set R^_i after 
removing these rays. Then 

rER° 

Furthermore, since i?? c and \ -R^'l < 3, 

Y.-l - = Y Pr-^"'^ - E (Pr-' - Pr-^'-'^) 

Further upper bounding by dropping the second sum shows that the resulting first term is still smaller 
than the bound (3.33). 
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As a result, we consider the larger bound (3.33) of these two cases and compute 

k —2k 
1=1 

Applying Theorem 3.4 completes the proof. □ 
Remark 3.3. Expanding the r.h.s. of (3.29) around in terms of the variable shows 

Pr [|X-7V^| >5] <2exp(^-— j for d^oo. (3.34) 

This indicates appropriate choices k = k{d) for large but finite problem sizes d occurring in applications, 
so as to bound the deviation of from its expected value. As a result, for such choices of k, our analysis 
based on expected values of the key system parameters will hold in applications with high probability. 

3.4. Critical Sparsity Values and Recovery. We derived the expected number NR{k) of nonzero measure- 
ments rrired (2.5) induced by random fc-sparse vectors x G ^ and the corresponding expected number 
NR(k) of non redundant cells rired- The tail bound. Prop. 3.5, guarantees that the dimensions of reduced 
systems concentrate around the derived expected values, explaining the threshold effects of unique recovery 
from sparse tomographic measurements. 

We now introduce some further notations and discuss the implication of Section 2.2 on exact recovery of 
X G Let Nji(k) and Nc{k) be the expected dimensions of the reduced system induced by a random 

fc-sparse nonnegative vector as detailed in Corollary 3.3. Let 5 = and denote by ks^ kcrit and ki/s the 
on d dependent sparsity values which solve the equations 

NR{ks) = 6iNc{ks), (3.35) 

NR(kcr^t) = Nc{kcr^t). (3.36) 

Nnikopt) = iNcikopt). (3.37) 

Nnikys) = ^Ncikys)- (3.38) 

In what follows, the phrase with high probability refers to values of the sparsity parameter k for which 
random supports | supp(6)| concentrate around the crucial expected value Nr according to Prop. 3.5, thus 
yielding a desired threshold effect. 

Proposition 3.6. The system Ax — h, with measurement matrix A, admits unique recovery of k-sparse 
non-negative vectors x with high probability, if 

k<^=:h. (3.39) 
1 + c) 

For perturbed systems we have. 

Proposition 3.7. The system Ax — b, with perturbed measurement matrix A, admits unique recovery of 
k-sparse non-negative vectors x with high probability, ifk satisfies condition k < kcrit- 

In case Conjecture 2.5 holds, uniqueness of x G M^j^^ + would be guaranteed if A: < ki/§. Finally, 
we comment on the maximal sparsity threshold kopt, in case reduced systems would follow a symmetric 
distribution with respect to the origin and columns would be in general position. 
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Figure 3 . Imaging setup with 4 cameras corresponding to the image planes, shown as two 
pairs in the left and center panel, respectively. Right panel: Cell centers projected onto the 
first image plane are shown as dots for the case d = 5. The cube = [0, rf]^ is discretized 
into cells and projected along 4 • d{2d — 1) rays. 



4. 4 Cameras - Left Degree equals 4 

We consider the imaging set-up depicted by Figure 3 and conduct a probabilistic analysis of its recovery 
properties, analogous to 3.2. This scenario is straightforward to realize and should also be particularly 
relevant to practical applications. 

4.1. Imaging Geometry. Each coordinate of the unit cube = [0, d]^ is discretized into the intervals 
{0, 1, 2, . . . , d}, resulting in d^ voxels with coordinates 

C={c= I) - ^(1, 1, 1) : i,j, le[d]}. (4.1) 

There are 4 sets of parallel projection rays corresponding to the normals of the image planes depicted in 
Fig. 3, 

= ^(-1, 0, 1), = i=(i, 0, 1), = -^(0, -1, 1), = -^(0, 1, 1). (4.2) 

We denote the set of projection rays and its partition corresponding to the 4 directions 

R^uURi. (4.3) 

Each set Ri contains {2d — 1) • d projection rays whose measurements yields a projection image with {2d — 
1) X d pixels. We index and denote the pixels by (5, t), and the projection rays through these pixels by 

zG {1,2,3,4}. 

For each cell c G C indexed by z, j, / G [d] according to (4.1), we represent the corresponding pixels after a 
suitable transformation by 

(51, ti) = {i + l-l-dj), si G [l-d,d-l], ti G [d], (4.4a) 

(52,^2) = S2 G [l-d,d-l], t2 G [d], (4.4b) 

(53,^3) = (i,j + /-l-d), 53 G [d], ts G [l-d,d-l], (4.4c) 

(54,^4) = {ij - /), 54 G [d], t4 G [1 - d, d- 1]. (4.4d) 
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The cardinahties of the projection rays, i.e. the number of cells covered by each projection ray, are 

aG{l,2}: \rl^\ = d - \s\, t G [d], (4.5a) 



6g{3,4}: |r^J = d-|t|, se [d], (4.5b) 



We observe the symmetries 



and define 



Ks,t\ = Kt\, Ks\ = \rU (4-6) 



^.1 Ksil (4-7) 



because [r^ j| does not vary with t. Summing up the cells covered by all rays along the first direction, for 
example, we obtain by (4.5a), (4.6) and (4.7), 

r^eRi te[d]s=l-d s=l-d s=l 



We set 



R{ki,k2) := 1 



Rik,, k,, k,) (l - Mj±y_LM^)\ (4.8) 



d3 ' 



We conduct next for this setup the analysis analogous to Section 3.2, in order to compute the expected 
size of the reduced system (2.5) for random fc-sparse vectors x. 

4.2. Dimensions of Reduced Systems. We first compute the expected number of measurements rrired (2.5) 
as a function of the sparsity parameter k. 

Lemma 4.1. The expected number ringed — of non-zero measurements is 

Nr = Nnik) = E[| supp(6)|] = \R\ - < = Ad{2d - 1) - N% (4.9a) 

Proof. Taking into account symmetry, we have 

r^R r^R r^eRi ' ' 

Applying (4.5a) yields the assertion. □ 

Proposition 4.2. The expected size rired — of subset of cells that support random subsets R^, C of 
observed non-zero measurements, is 

Nc = Nc{k) ^d^ -N}. + NI-NI + N^ (4.10) 
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where 

d-i 



^ S=l ^ 



]^l^2d^ R(l^i-\-d,i-l)^\ ^ R{l-i,l-j), 

i.leld] ij^l^ld] (4.11) 

N^ = 2 {R{l + i-l-d,l-i,l-j) + R{l-i,l-j,l + j-l-d)), 

i,j,le[d] 

R{l + i- l-d,l-i,l + j - l-d,l-j), 

ij,le[n] 

end the functions R are given by (4.8). 

Proof. We consider for each cell c G C the quadruple of projection rays (r^, r^y^r^^ Tq) meeting in this cell, 
and the corresponding partition (4.3) of projection rays. Cell c is contained in the set Ch (2.4) supporting i?5 
if not any ray of the corresponding quadruple returns a zero measurement. Thus, 



TVc = e[5^(1 - X,i)(l - X,2)(l - X,3)(l - X,. 

4 



cGC i=l l<i<j<4 l<i<j<K4 



l<^<j<K4 



d3 ; V 



(4.12) 



We consider each term in turn. 

(i) As for the first term, we obviously have \C\ = d^. 

(ii) Concerning the second term, taking symmetry into account we compute, 

EE(-^r-E(-tpr^.EE(-^^' 

c^C i=l c^C r^eRicer^ 

Since = for all c G r^, we obtain using (4.5), 

ceC i=l r^&Ri si=l-d 

1 Nfc '^"^ 

= M{ d{ 1 



(iii) Concerning the third term of (4.12), we consider first the contribution of the pair of directions (z, j) = 
(1,2). Replacing cell-indices of projection rays by pixel indices according to (4.4), we have using (4.5) 
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and (4.7), 

-dj2[i- h±hiA±hA^y (4.13) 

i,le[d] 

= d R{i + l-l-d,i-l), 

i.leld] 

where the factor d appears because the summand does not depend on j by (4.5a), and R is defined by 
(4.8). For the pair of directions (3, 4), we get 

(l - - E (l - ^ -1^^ ) ' (4.14) 

which equals (4.13) due to the symmetry (4.6). 

Next, we consider the pair of directions (1,3). Taking into account the symmetry (4.6) and using 
(4.8), we obtain 

^}Ur^\\k / \^ii+i-i-d,i)\^\^fu+l-i-d)\~'^\^ 



d^ 

^ ^ 1^-^ h+i-i-dl + \rj+i-i-d\ - ^ 1^-^ _ In-il + In-jl - 



d^ 



In the same way it can be shown that the remaining pairs of directions (1,4), (2, 3), (2, 4) each con- 
tributes the last expression, 
(iv) Concerning the fourth term of (4.10), we get for the triple of directions (1, 2, 3) the contribution 



^ j^^ \rl^r^c^r^\ y ^ ^ \rt+i-i-d\ + \r^-l\ + \rj+i-i-d\ - 

= i?(z + /-l-rf,/-z,j + /-l-d)), 

ij,le[n] 

and likewise for the remaining triples 

(1,2,4): Yl Rii + l-l-dJ-iJ-j), 

ij,le[n] 

(1,3,4): J2 Rii + l-l-dJ + l-l-dJ-j), 

ij,le[n] 

(2,3,4): J2 R{l-iJ + l-l-d,l-j). 

Evidently, the first and last pair of expressions are equal, respectively, 
(v) Finally, the expression for the last term of (4.12) is immediate. 



(4.16) 



(4.17) 



□ 
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Figure 4. Success and failure empirical phase transitions for the 2D, 3-cameras case, 
hexagonal area, from Section 3, Fig. 2, left. Reduced unperturbed (top left) and perturbed 
(top right) matrices are overdetermined and of full rank with hight probability if the cor- 
responding sparsity level is below ks (unperturbed case) or kcrit (perturbed case). The 
<-marked curve depicts fc^/d^ (3.35), and the D-marked curve, kcrit /d? from (3.35). Prob- 
ability of uniqueness in [0, 1]^ of a A: = pdP sparse binary vector for unperturbed (bottom 
left) and perturbed matrices (bottom right). This probability is high below ks/d? and de- 
creases slowly, for the unperturbed case (bottom left). In the perturbed case the empirical 
probability of uniqueness exhibits a sharp transition accurately described by the [>-marked 
curve kiis/d? from (3.38), which for the 3 camera case lies under the o-marked curve kopt 
from (3.37). 



We conclude this section by stressing that critical sparsity values ks (3.35), kcrit (3.36), kopt (3.37), kij^) 
(3.38), can be worked out based on the just derived values Nr and Nq- A tail bound may be derived 
analogously to Prop. 3.5. We omit this redundant detail due to space constraints. 

5. Numerical Experiments and Discussion 

In this section we relate the previously derived bounds on the required sparsity that guarantee unique 
nonnegative or binary fc-sparse solutions to numerical experiments. In analogy to [6] we assess the so 
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Figure 5. Left: The analytically derived curves from Section 3 correctly follow the con- 
tour lines of the average fraction of reduced systems determined empirically as a function of 
d and relative sparsity k/d. From bottom to top: ks/d (<-marked curve), kcrit / d (D-marked 
curve), kopt/d (o-marked curve) and ki/^/d (<]-marked curve). Right: These curves are 
plotted again (light gray) for a wider range and compared to the analogous curves for the 
geometry in Fig. 2, right, square area, 2 orthogonal, one diagonal projecting direction. 

called phase transition p as a function of d, which is reciprocally proportional to the undersampling ratio 
^ G (0, 1). We vary d, build a specific matrix projection matrix A along with its perturbed version A and 
consider the sparsity as a fraction of d in 2D or in 3D, respectively, thus k = pd^~^ , with p G (0, 4) and 
D G {2, 3}. This phase transition p{d) indicates the necessary relative sparsity to recover a fc-sparse solution 
with overwhelming probability. More precisely, if ||x||o < p{d)-d^~^, then with overwhelming probability a 
random fc-sparse nonnegative (or binary) vector x* is the unique solution in := {x : Ax = Ax*, x > 0} 
or Jo,i •= Ax = Ax*,^ G [0, 1]^}, respectively. Uniqueness can be "verified" by minimizing and 
maximizing the same objective f^x over or Tq^i, respectively. If the minimizers coincide for several 
random vectors / we claim uniqueness. The resulting linear programs we solved by a standard LP solver 
^ As shown e.g. in Fig. 9 and confirmed by all our numerical experiments the threshold for a unique 
nonnegative solution and a unique 0/1-bounded solution are quite close, especially for high values of d. 

We note that A has the same sparsity structure as A, but random entries drawn from the standard uniform 
distribution on the open interval (0.9, 1.1). 

Then for p G [0, 1] a pd^^^-sparse nonnegative or binary vector was generated to compute the right 
hand side measurement vector and for each (d, p) -point 50 random problem instances were generated. A 
threshold-effect is clearly visible in all figures exhibiting parameter regions where the probability of exact 
reconstruction is close to one and it is much stronger for the perturbed systems. The results are in excellent 
agreement with the derived analytical thresholds. We refer to the figure captions for detailed explanations and 
stress that a threshold-effect is clearly visible in all figures exhibiting parameter regions where the probability 
of exact reconstruction is close to one. 

5.1. 2D: Hexagonal Volume and 3 Cameras. In this 2D, 3-cameras case, we generated A according to 
the geometry described in Section 3, Fig. 2, left. We considered d G {51, 71, 91, ... , 251} and varied the 
sparsity k = pd, by varying p in (0, 0.5) with constant stepsize 0.01. The obtained empirical phase transitions 
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Figure 6. Empirical relative critical curves in for 3 to 8 cameras in 2D. Top left: ks/d 
(>-marked curve). Top right: kcrit/d (D-marked curve). Bottom left: kopt/d (o-marked 
curve). Bottom right: fci/^/rf (<-marked curve). Perfect recovery is possible for sparsity 
levels below ki/s for perturbed systems. 

are depicted in 4. Fig. 7 additionally shows how the analytically determined critical curves compare to the 
critical curves obtained empirically for the geometry in Fig. 2, right. The (4(i — 1) x (P projection matrix 
corresponding to a square area and three projecting directions has similar reconstruction properties and the 
critical curves slightly change by factor ^. 

5.2. 2D: Square Volume and 3 to 8 Cameras. In our theoretical analysis in the previous sections we 
derived the expected number of nonzero rows Nji{k) induced by the fc-sparse vector along with the num- 
ber Nc{k) of "active" cells which cannot be empty. This can be done also empirically, compare Fig. 5, 
left. We have done this in 2D up to 8 projecting directions and obtained empirically the critical curves 
ks, kcriu Kpt and ki/s- To generate the curves we varied k/d G (0,4) by a constant stepsize 0.01 and 
d G {50, 100, • • • , 1050} and generated for each point (fc/d, d) 500 problem instances. Further we de- 
termined the contour lines of Nji{k/d, d)/Nc{k/d, d) corresponding to the levels {5^, 1, 0.5, ^7^}. The 
relative sparsity curves ks/d, kcrit/d, kopt/d and ki/^/d are plotted in Fig. 6 and accurately follow the 
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d=200 




Figure 7. Left: kopt/d (o-marked curve) along with ki/s/d (<]-marked curve). For 3 cam- 
eras ki/s/d lies below kopt/d, but starting with 5 cameras ki/s significantly outperforms 
kopt' This shows the fundamental difference between the considered 0/1 -matrices and ran- 
dom matrices underlying a symmetrical distribution with respect to the origin. For ran- 
dom matrices recovery of fc-sparse positive or binary vectors with sparsity levels beyond 
kopt would be impossible. Right: For d = 200 the 6 curves depict the average ratio of 
^red{k) / riredik) as a function of sparsity k for 3 to 8 cameras from bottom to top. 



# 


m 


n 


projection angles 


3rd camera 


4d- 1 




0°,90°,45° 


4th camera 


6d-2 




0°,90°,t45° 


5th camera 


7d+ [fj -2 




0°,90°,T45°,arctan(2) 


6th camera 


8d + 2[fj -2 




0°,90°,^45°,^arctan(2) 


7th camera 


9d + 3[fj -2 




0°, 90°, ^45°, ^ arctan(2), arctan(0.5) 


8th camera 


I0d + 4[|J -2 




0°, 90°, ^45°, ^ arctan(2), ^ arctan(0.5) 



Table 1 . Dimensions of full projection matrices. 



empirical recovery thresholds, as shown e.g. in Fig. 8 for 6 cameras. Fig. 7 shows Nji^k, 200) /Nc{k, 200) 
for varying number of cameras. The projection angles we chosen such that the intersection with all cells 
is constant, yielding binary projection matrices after scaling. Each camera resolution differs with different 
angle. We summarize the used parameters in Table 1. 

5.3. 3D: 4 Cameras. In 3D we consider the matrix from Section 4, Fig. 3, vary d G {25, 26, ... , 55} and 
k = pd'^ by varying p G (0,4) with stepsize 0.01. For larger values of d the empirical thresholds follow 
accurately the estimated curves. Note that for = 55, A is a 48180 x 166375 matrix. 

6. Conclusion 

The new measurement paradigm of compressed sensing seeks to capture the "essential" aspects of a high- 
dimensional but sparse object using as few measurements as possible by randomization. Provided that the 
measurements satisfy certain properties nonnegative sparse signals can be reconstructed exactly from a sur- 
prisingly small number of samples. Moreover, there exist precise thresholds on sparsity such that for any 
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Figure 8 . Success and failure empirical phase transitions for the 2D, 6 cameras case. Prob- 
ability of uniqueness in [0, 1]^ of a = pd? sparse binary vector for unperturbed (left) and 
perturbed matrices (right), along with ks/d^ (t>-marked curve), kcHt/d^ (D-marked curve), 
kopt/d^ (o-marked curve), and fci/^/rf^ (<-marked curve) from bottom to top. The dashed 
line depicts k/dP, with k solving mred{k) = jTiredik), which accurately follows the border 
of the highly success area/t^r all considered number of cameras, 3, 4 ... 8. Recovery is pos- 
sible beyond kopt, accurately described by ki/^. In the 6 cameras case there is no evident 
performance boost for perturbed systems, since the columns of reduced systems are most 
likely to be in general position for both perturbed and unperturbed systems. However, in the 
perturbed case, recovery is more stable. 

nonnegative solution that is sparser than the threshold is also the unique nonnegative solution of the under- 
lying linear system. Tomographic projection matrices do not satisfy the conditions which allow applying 
these results on the image reconstruction problem even if the sought solution is very sparse. However, we 
analytically showed in the present work that there are thresholds on sparsity depending on the numbers of 
measurements, below which uniqueness is guaranteed and recovery will succeed and above which it fails 
with high probability. When recovery succeeds it yields perfect reconstructions, without any ghost-particles. 
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Figure 10. Success and failure empirical phase transitions for the 3D, 4 cameras case, 
from Section 4, Fig. 3. Probability of uniqueness in [0, 1]^ of a A: = pd? sparse binary 
vector for unperturbed (top left) and perturbed matrices (top right). The t>-marked curve 
depicts ks/(f (3.35), the D-marked curve kcHt/d^ (3.36), the o-marked curve kopt (3.37) 
and the <-marked curve kij^/dP (3.38). In case of the perturbed matrix A exact recovery 
is possible beyond kopt/d?- Moreover k^i^l^ follows most accurately the empirical phase 
transition for perturbed systems for high values of d. The empirical probability that the 
reduced unperturbed matrices are overdetermined and of full rank (bottom figure), exhibits 
a threshold in between the estimated relative critical sparsity level k^ and kcrit- This explains 
the performance boost for perturbed systems. 



